Когда признаков много, можно представить все объекты как облако точек в многомерном пространстве.
Migration by Don McCullough on Flickr
black shadows for a white horses / les negres ombres dels cavalls blancs by Ferran Jordà on Flickr
Метод визуализации отношений между объектами в пространстве с небольшим числом измерений (обычно 2).
Исходные данные — матрица расстояний между объектами в многомерном пространстве.
nMDS ординация старается сохранить отношения между объектами.
Расстояния по автородогам не совсем “евклидовы”: например, из Кронштадта до Санкт-Петербурга нельзя добраться по прямой.
Карта
Проекция при помощи анализа главных координат (PCoA)
Карта
Ординация nMDS
Ординация сохраняет ранг расстояний между объектами (похожие располагаются близко, непохожие — далеко; если объект А похож на B, больше чем на C, то и на ординации он, скорее всего, окажется ближе к B, чем к C).
На ординации имеет смысл только взаиморасположение объектов. Облако точек в осях MDS можно вращать, перемещать, зеркально отражать. Суть ординации от этого не изменится.
Значения координат точек в ординации лишены смысла (их вообще можно не приводить на итоговой ординации).
Взаиморасположение точек на плоскости подобно взаиморасположению точек в многомерном пространстве признаков, но значения не полностью совпадают.
Мы можем построить график зависимости расстояний между точками на ординации и соответствующих им значений коэффициентов различия в исходном пространстве.
У хорошей ординации nMDS ранги расстояний должны совпадать с рангами коэффициентов различия в исходном многомерном пространстве.
Диаграмма Шепарда показывает соответствие расстояний на ординации и коэффициентов различия в исходном многомерном пространстве.
В идеале, большим расстояниям в исходном пространстве должны соответствовать большие расстояния на ординации, а меньшим — меньшие (т.е. должны сохраняться ранги расстояний).
Красная линия — предсказанные монотонной регрессией значения расстояний на ординации — монотонно возрастает.
Монотонная регрессия минимизирует сумму квадратов отклонений значений координат на ординации от значений в исходном многомерном пространстве, сохраняя монотонные отношения между ними.
Алгоритм подбора монотонной регрессии называется Pool-Adjacent-Violator-Algorithm (PAVA)
Стресс (Stress) показывает, насколько соответствуют друг другу взаиморасположение точек на плоскости ординации и в исходном многомерном пространстве признаков.
Стресс — это корень из суммы квадратов остатков монотонной регрессии, отнесенной к сумме квадратов всех коэффициентов различия в исходной матрице.
\[ Stress(1) = \sqrt{\frac{\sum_{h < i}{(d_{h,i} - \hat d_{h,i})^2}}{\sum_{h < i}{d_{{h,i}}^2}}}\]
Хорошо, когда \(Stress(1) < 0.2\).
Данные можно скачать с сайта Pangaea
dat <- read.delim("data/Krapivin_2017_Medulis-symb.tab", skip = 36)
str(dat)
## 'data.frame': 548 obs. of 23 variables: ## $ Event : Factor w/ 3 levels "Gorelyi_Isl",..: 3 3 3 3 3 3 3 3 3 3 ... ## $ Date.Time : Factor w/ 1 level "2013-07": 1 1 1 1 1 1 1 1 1 1 ... ## $ Ord.No : int 1 2 3 4 5 6 7 8 9 10 ... ## $ Area : Factor w/ 3 levels "B. Gorelyi isl.",..: 2 2 2 2 2 2 2 2 2 2 ... ## $ Site : int 1 1 1 1 1 1 1 1 1 1 ... ## $ Zone : Factor w/ 3 levels "down","middle",..: 1 1 1 1 1 1 1 1 1 1 ... ## $ Comment : int 1 1 1 1 1 1 1 1 1 1 ... ## $ M..edulis.L.shell..mm.: num 31.2 41.1 46.4 40 34.5 47.6 34.3 49.6 49.8 30 ... ## $ M..edulis.age..a. : num 1 2 2 3 4 3 2 4 3 1 ... ## $ Urastoma.spp..... : int 26 16 19 11 14 16 10 4 11 0 ... ## $ Renicola.spp..... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Himasthla.spp..... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Metacercaria.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Gymnophallus.spp..... : int 0 0 0 0 0 0 0 0 4 0 ... ## $ Bucephalidae.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Alg.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Nematoda.... : int 1 1 0 0 0 0 0 0 0 0 ... ## $ Microsetella.... : int 0 0 0 0 0 2 0 0 0 0 ... ## $ Copepoda.... : int 0 0 0 0 0 1 0 0 0 0 ... ## $ Chironomidae.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Halacaridae.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Jaera.spp..... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Ostrac.... : int 0 0 0 0 0 0 0 0 0 0 ...
colnames(dat) # Было
## [1] "Event" "Date.Time" "Ord.No" ## [4] "Area" "Site" "Zone" ## [7] "Comment" "M..edulis.L.shell..mm." "M..edulis.age..a." ## [10] "Urastoma.spp....." "Renicola.spp....." "Himasthla.spp....." ## [13] "Metacercaria...." "Gymnophallus.spp....." "Bucephalidae...." ## [16] "Alg...." "Nematoda...." "Microsetella...." ## [19] "Copepoda...." "Chironomidae...." "Halacaridae...." ## [22] "Jaera.spp....." "Ostrac...."
colnames(dat) <- gsub("[.]", replacement = "", colnames(dat))
dat <- rename(dat, L = "MedulisLshellmm", Age = "Medulisagea")
colnames(dat) # Стало
## [1] "Event" "DateTime" "OrdNo" "Area" ## [5] "Site" "Zone" "Comment" "L" ## [9] "Age" "Urastomaspp" "Renicolaspp" "Himasthlaspp" ## [13] "Metacercaria" "Gymnophallusspp" "Bucephalidae" "Alg" ## [17] "Nematoda" "Microsetella" "Copepoda" "Chironomidae" ## [21] "Halacaridae" "Jaeraspp" "Ostrac"
# Делаем сайт фактором
dat$Site <- factor(dat$Site, levels = c(2, 3, 1), labels = c("G","L","S"))
# Сливаем редкие виды в одну категорию
f_remove <- c("Nematoda", "Microsetella", "Copepoda",
"Chironomidae", "Halacaridae", "Jaeraspp",
"Ostrac")
dat$Other <- rowSums(dat[, f_remove])
# Суммарная численность симбионтов
f_sp <- c("Urastomaspp", "Renicolaspp", "Himasthlaspp",
"Metacercaria", "Gymnophallusspp", "Alg",
"Other")
dat$Total <- rowSums(dat[, f_sp])
# Данные для анализа
# Только мидии с симбионтами и возрастом от 3 до 8 лет
dat <- dat[dat$Total != 0 & dat$Age %in% 3:8, ]
spec <- dat[, f_sp] # виды-симбионты
env <- dat[, c("Zone", "Site", "L", "Age")] # свойства мидий-хозяев
Постройте ординацию мидий по обилию разных видов-симбионтов с использованием коэффициента различия Брея-Куртиса. Следите, чтобы алгоритму удалось найти финальное решение. Если необходимо, настройте вызов metaMDS()
Вычислите стресс для получившейся ординации.
Нарисуйте простейший график.
ord_mus <-
ord_mus <- metaMDS(spec, dist = "bray")
По-умолчанию metaMDS() пытается подобрать оптимальную трансформацию, подходящую для анализа сообществ. В данном случае — извлекает корень и делает двойную висконсинскую стандартизацию.
Square root transformation Wisconsin double standardization Run 0 stress 0.1319177
Такая модель не сходится.
*** No convergence -- monoMDS stopping criteria:
19: stress ratio > sratmax
1: scale factor of the gradient < sfgrmin
После отключения автоматического подбора трансформации данных удается найти решение.
ord_mus <- metaMDS(spec, dist = "bray", autotransform = FALSE)
## Run 0 stress 0.141 ## Run 1 stress 0.145 ## Run 2 stress 0.141 ## ... Procrustes: rmse 0.00388 max resid 0.0344 ## Run 3 stress 0.145 ## Run 4 stress 0.146 ## Run 5 stress 0.144 ## Run 6 stress 0.154 ## Run 7 stress 0.154 ## Run 8 stress 0.141 ## ... Procrustes: rmse 0.00221 max resid 0.0231 ## Run 9 stress 0.146 ## Run 10 stress 0.141 ## ... Procrustes: rmse 0.00273 max resid 0.0447 ## Run 11 stress 0.141 ## ... Procrustes: rmse 0.0034 max resid 0.0459 ## Run 12 stress 0.141 ## ... Procrustes: rmse 0.00311 max resid 0.034 ## Run 13 stress 0.144 ## Run 14 stress 0.141 ## ... Procrustes: rmse 0.00342 max resid 0.0419 ## Run 15 stress 0.141 ## ... Procrustes: rmse 0.00299 max resid 0.038 ## Run 16 stress 0.141 ## ... Procrustes: rmse 0.00318 max resid 0.0343 ## Run 17 stress 0.144 ## Run 18 stress 0.15 ## Run 19 stress 0.141 ## ... Procrustes: rmse 0.00214 max resid 0.0217 ## Run 20 stress 0.141 ## ... Procrustes: rmse 0.00236 max resid 0.0329 ## *** No convergence -- monoMDS stopping criteria: ## 20: stress ratio > sratmax
Хорошее решение, но график нужно украсить.
ord_mus$stress
## [1] 0.141
ordiplot(ord_mus)
# Палитры
pal_col <- c("red", "green", "steelblue")
pal_sh <- c(1, 2, 0)
# Украшенный график
ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])
ordiplot(ord_mus, type = "n") points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site]) text(ord_mus, display = "spec", cex = 0.9, col = "grey20")
Используя данные, приведенные ниже, постройте график ординации при помощи ggplot2. Покажите точками — мидий, текстом — центроиды для видов-симбионтов, цветом — зону литорали, формой — сайт, размером маркеров — размер мидий.
# Координаты точек (мидий) ord_mus_pt <- data.frame(env, scores(ord_mus, display = "sites")) head(ord_mus_pt, 2)
# Координаты центроидов переменных (видов-симбионтов) ord_mus_sp <- data.frame(scores(ord_mus, display = "species")) ord_mus_sp$Species <- rownames(ord_mus_sp) head(ord_mus_sp, 2)
gg_ord_mus <- ggplot() +
geom_point(data = ord_mus_pt,
aes(x = NMDS1, y = NMDS2, colour = Zone,
shape = Site, size = L), alpha = 0.5)
gg_ord_mus
gg_ord_mus_sp <- gg_ord_mus +
geom_text(data = ord_mus_sp,
aes(x = NMDS1, y = NMDS2,
label = Species))
gg_ord_mus_sp
envfit()Функция envfit() строит регрессионную модель зависимости переменной среды от координат точек на плоскости ординации.
Модель вот такого вида:
\(y_i = \beta_1 \text{NMDS1}_i + \beta_2 \text{NMDS2}_i + \varepsilon_i`\)
Значимость зависимости оценивается при помощи пермутационного теста (это позволяет работать в ситуации, когда неполностью выполнены условия применимости).
Такой вызов функции позволит подобрать сразу несколько моделей и изобразить их на графике ординации:
ef <- envfit(ord_mus, env[, c("Zone", "Site", "L", "Age")])
envfit() для непрерывных переменныхef$vectors
## NMDS1 NMDS2 r2 Pr(>r) ## L -1.000 -0.030 0.19 0.001 *** ## Age 0.790 -0.613 0.04 0.001 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## Permutation: free ## Number of permutations: 999
MDS1 и MDS2 содержат информацию о направлении вектора непрерывной переменной. Длина вектора пропорциональна корню из коэффициента детерминации для модели.r2 — \(R^2\) из регресcионного анализа.Pr(>r) — пермутационная оценка статистической значимости.envfit() для дискретных переменныхef$factors
## Centroids: ## NMDS1 NMDS2 ## Zonedown -0.65 -0.09 ## Zonemiddle -0.14 -0.05 ## Zoneup 0.83 0.14 ## SiteG 0.19 0.17 ## SiteL 0.23 -0.41 ## SiteS -0.42 0.22 ## ## Goodness of fit: ## r2 Pr(>r) ## Zone 0.37 0.001 *** ## Site 0.16 0.001 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## Permutation: free ## Number of permutations: 999
MDS1 и MDS2 содержат координаты центроидов для дискретных факторов (т.е. положение центра облака точек, относящихся к категории фактора).r2 – \(R^2\) из регресcионного анализа.Pr(>r) — пермутационная оценка статистической значимости.envfit()ordiplot(ord_mus, type = "n") points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site]) plot(ef)
# install.packages("devtools")
# devtools::install_github("gavinsimpson/ggvegan")
library(ggvegan)
ord_mus_ef <- fortify(ef)
ord_mus_ef
envfit()gg_ord_mus +
geom_segment(data = ord_mus_ef[ord_mus_ef$Type == "Vector", ],
aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm"))) +
geom_text(data = ord_mus_ef[ord_mus_ef$Type == "Vector", ],
aes(x = NMDS1, y = NMDS2, label = Label, hjust = 1.1, vjust = 1)) +
geom_text(data = ord_mus_ef[ord_mus_ef$Type == "Centroid", ],
aes(x = NMDS1, y = NMDS2, label = Label, hjust = 1.1, vjust = 1))
envfit()envfit() основан на предположении о линейном характере связи между внешней переменной и координатами на ординации. Это не всегда так.ordisurf()ordisurf() применяется, когда зависимость нелинейна.
В основе работы функции — общая аддитивная модель (General Additive Model, GAM).
GAM позволяют подобрать нелинейные зависимости, даже если заранее неизвестна их форма.
ordisurf() использует пакет mgcv для подбора GAM на основе метода тонких пластин (thin plate spline).
ordisurf()Изолинии на графике ординации показывают предсказанное значение зависимой переменной для разных участков ординации.
par(mfrow = c(1, 2)) os_L <- ordisurf(ord_mus, env$L, method = "REML") os_Age <- ordisurf(ord_mus, env$Age, method = "REML") par(mfrow = c(1, 1))
ordisurf()Длина мидий линейно меняется на ординации (эффективное число степеней свободы для сплайна \(edf=1.96\)). Зависимость довольно хорошо выражена (объясненная девианса 18.7%).
summary(os_L)
## ## Family: gaussian ## Link function: identity ## ## Formula: ## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE) ## ## Parametric coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 43.410 0.484 89.7 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df F p-value ## s(x1,x2) 1.96 9 9.83 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## R-sq.(adj) = 0.183 Deviance explained = 18.7% ## -REML = 1457.6 Scale est. = 92.439 n = 395
ordisurf()Возраст мидий меняется нелинейно вдоль ординационных осей (эффективное число степеней свободы для сплайна \(edf=4.99\)). Эта зависимость довольно слабая (объясненная девианса 8.97%).
summary(os_Age)
## ## Family: gaussian ## Link function: identity ## ## Formula: ## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE) ## ## Parametric coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.8861 0.0713 68.5 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df F p-value ## s(x1,x2) 5.06 9 3.76 0.00000051 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## R-sq.(adj) = 0.0791 Deviance explained = 9.09% ## -REML = 704.9 Scale est. = 2.0097 n = 395
fortify.ordisurf <- function(model) {
# Fortifies an object of class ordisurf to produce
# a dataframe of contour coordinates
# for subsequent plotting in ggplot.
xy <- expand.grid(x = model$grid$x, y = model$grid$y)
xyz <- cbind(xy, c(model$grid$z))
names(xyz) <- c("x", "y", "z")
return(na.omit(xyz))
}
ord_mus_os <- fortify.ordisurf(os_Age)
head(ord_mus_os, 4)
ordisurf()ggplot(data = ord_mus_os, aes(x = x, y = y, z = z)) + stat_contour(aes(colour = ..level..)) + labs(x = "NMDS1", y = "NMDS2", colour = "Age")
ordisurf()ordisurf() важно, чтобы в ординации участвовало много объектов и они располагались на плоскости ординации более-менее равномерно. Пустоты и отскоки снижают качество предсказаний.
envfit() и ordisurf() обращайте внимание на форму связи.Во всех примерах необходимо:
Источники данных
ade4).vegan).